### Need to create some plots to examine microCLIM
### Looking for how it varies with respect to BIOCLIM
### And how it varies with respect to topographic features
### Start by assembling a dataframe of microCLIM, BIOCLIM, and topographic features for all grid cells in the WTWHA

### Load library

library('SDMTools')
library('akima')

### Establish directories

bc.dir = '/home1/99/jc152199/MicroclimateStatisticalDownscale/250mASCII/BC6_monthly/'
mc.dir = '/home1/99/jc152199/MicroclimateStatisticalDownscale/250mASCII/NEWmicroCLIM_76-05/'
topo.dir = '/home1/99/jc152199/MicroclimateStatisticalDownscale/250mASCII/STATIC/'
#out.dir = '/home1/99/jc152199/microCLIM/images/'
out.dir = '/home1/99/jc152199/microCLIM/testimages/'

### List ASCII files of interest in the above directories

bc.files = list.files(bc.dir, full.names=T, pattern='.asc.gz')[1:4]
mc.files = list.files(mc.dir, full.names=T, pattern='.asc.gz')
topo.files = list.files(topo.dir, full.names=T)[c(1,3:length(list.files(topo.dir)))]
all.files = c(bc.files,mc.files,topo.files)

### List names for the above files

bc.names = gsub('.asc.gz','',list.files(bc.dir, pattern='.asc.gz')[1:4])
mc.names = gsub('_','',gsub('.asc.gz','',list.files(mc.dir, pattern='.asc.gz')))
topo.names = gsub('_WTplusbuffer_LatLong_WGS1984_250mres.asc','',list.files(topo.dir)[c(1,3:length(list.files(topo.dir)))])
all.names = c(bc.names,mc.names,topo.names)

### Loop through and read in these files into a single dataframe

out.data = NULL ### Null object to bind data onto

for (i in 1:length(all.files))

	{
	
	if(i<9) {t.asc = as.vector(read.asc.gz(paste(all.files[i],sep='')))}
	if(i>8) {t.asc = as.vector(read.asc(paste(all.files[i],sep='')))}
	
	out.data = cbind(out.data,t.asc)

	cat('\n',all.names[i],' - Processing Completed\n',sep='')
	
	}
	
# Close loop, dataframe assembled

### Convert out.data to a dataframe

out.data = data.frame(out.data)

### Change names of out.data

names(out.data) = all.names

### Remove rows where DEM is 'NA'

out.data = out.data[which(is.na(out.data$dem)==F),]

### Write out out.data

write.csv(out.data, file=paste(out.dir,'ASCII_dataframe.csv',sep=''),row.names=F)

### Read in out.data

out.data = read.csv('/home1/99/jc152199/microCLIM/images/ASCII_dataframe.csv',header=T)

### For faster plotting randomly sample 1/50th of out.data while preparing this script

out.data = out.data[c(sample(c(1:nrow(out.data)),nrow(out.data)/50,replace=FALSE)),]

### Start by plotting up BC's versus their respective MC's

for (i in 1:4)

	{
	
	png(paste(out.dir,'Surface0',substr(names(out.data)[i],4,4),'_Comparison.png',sep=''),units='in',height=6,width=6, res=300) # Open .png device driver
	
	plot(out.data[,i+4],out.data[,i],xlab='microCLIM',ylab='BIOCLIM',main = paste('Surface0',substr(names(out.data)[i],4,4),sep=''), cex=.5) # Plot values of interest
	
	lm1 = lm(out.data[,i]~out.data[,i+4]) # Perform a linear regression using microCLIM as the predictor and BIOCLIM as the predictant
	
	abline(a=0,b=1,col='red',lty=2,lwd=2) # Plot a 1:1 line
	
	abline(a=lm1$coefficients[1],b=lm1$coefficients[2],lty=1,lwd=2, col='blue') # Plot the line of the relationship between the surfaces
	
	legend('topleft',legend=c('1:1','Surface LM'), text.col = c('red','blue')) # Plot a legend
	
	dev.off()
	
	cat('\nSurface 0',substr(names(out.data)[i],4,4),' Plotting Complete\n',sep='')
	
	}
	
# Close plotting loop

### Now need to plot the BC/MC residuals in 2 and 3 dimensions
### Using the primary topographic drivers of the model as predictors of the residual values
### Calculate the difference between microCLIM and BIOCLIM, using microCLIM as the 'observed' and BIOCLIM as the 'predicted'

out.data$diff01 = out.data$microCLIM01-out.data$bc01
out.data$diff04 = out.data$microCLIM04-out.data$bc04
out.data$diff05 = out.data$microCLIM05-out.data$bc05
out.data$diff06 = out.data$microCLIM06-out.data$bc06

### Loop through and plot residuals against topographic variables

for (i in 9:15) # Loop through topographic variables first

	{
	
	png(paste(out.dir,names(out.data)[i],'_Vs_ClimateResiduals.png',sep=''),units='in',width=10, height=2.5,res=300) # Open .png driver
	
	par(mfrow=c(1,4)) # Configure plot space to a panel of 1 * 4
	
	for (j in 16:19)
	
		{
		
		plot(out.data[,i],out.data[,j],xlab=paste(names(out.data)[i],sep=''),ylab = paste(names(out.data)[j],sep=''),main = paste(names(out.data)[i],'_Vs_ClimateResiduals',sep=''), cex=.7)
		
		lm1 = lm(out.data[,j]~out.data[,i])
		
		abline(a=lm1$coefficients[1],b=lm1$coefficients[2],col='red',lty=2,lwd=2)
		
		cat('\n',names(out.data)[i],' versus ',names(out.data)[j],' comparison plotted\n',sep='')
		
		}
	
	dev.off()
	
	}
	
# Close plotting loop

### Now try a 3-D plot using DEM on the x-axis, FPC on the y-axis, and climate residual as the z-value

### Need remove NA's

out.data = na.omit(out.data)

### Create a 3 column matrix to plot the data in persp

zz = interp(out.data$dem,out.data$fpc,out.data$diff01,duplicate='strip')

### Create breaks for the colors
#min z value for diff01 = -4.5
#max z value for diff02 = +1.5

col.levs.below = seq(round(min(zz$z,na.rm=T)-1,0),0,.25)+.125
col.levs.above = seq(0,round(max(zz$z,na.rm=T)+1,0),.25)-.125
levs = c(col.levs.below,col.levs.above[3:length(col.levs.above)])

### Create a color pallete

cols.below = colorRampPalette(c('blue','green','white'))
cols.above = colorRampPalette(c('white','orange','red'))

col.codes.below = cols.below(length(col.levs.below)-1)
col.codes.above = cols.above(length(col.levs.above[3:length(col.levs.above)]))



### Now plot the object zz using the image command

png(paste(out.dir,'persptest.png',sep=''))
image(zz)
filled.contour(zz,levels=levs,col = c(col.codes.below,col.codes.above),xlab="dem",ylab="fpc")
dev.off()









	
	
